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ABSTRACT 


The  Trajectory  Analysis  Program  (TAP)  is  designed  specifically 
for  use  on  the  CDC  6000  series  machine  as  a  computational  tool 
to  generate  trajectories  by  different  mathematical  formulations 
and  to  compare  the  results.  In  particular,  the  reference  trajec¬ 
tory  is  the  result  of  numerical  integration  of  the  equations  of 
motion  in  a  Cowell  formulation  with  time  as  the  independent  vari¬ 
able;  whereas,  the  comparison  trajectory  selected  is  achieved  by 
an  analytic  (Pines,  Kyner,  or  Escobal)  or  a  numerically  integrated 
Encke  (classical  or  modified)  formulation.  A  comprehensive 
description  is  presented  of  the  algorithms  implemented  in  TAP 
with  emphasis  on  characteristics  such  as  logic  flow,  equations, 
numerical  techniques,  comparison  differencing,  and  plotting.  A 
basic  usage  guide  i3  included  with  complete  instructions  for  input 
data  preparation  and  explanations  of  output  tables. 
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NOMENCLATURE 


[A]  sum  of  constant  matrices  (evaluated  at  epoch)  whose  scalar 

coefficients  are  slowly  varying  functions  of  <p  and  the  rate 
parameters  (tjjTjV) 

[A]  time  derivative  of  the  matrix  [a] 

a  semimajor  axis 

a  mean  equatorial  radius  of  earth  or  semi-diameter  of 

e  the  central  mass 

a.  scalar  coefficients  of  rotation  matrix  [A] 

J 

CjyA/'W  drag  coefficient 

DEQ  numerical  integrator  used  for  trajectory  generation 

DP  double  precision 

ECI  earth- centered  inertial 

e  eccentricity 

F,  G  Keplerian  scalar  functions 

f  true  anomaly 

— LH  gravitational  acceleration  vector  due  to  noncentral  force 

field  of  earth  in  local  horizontal  system,  in  which 
coordinate  axes  are  directed  Up  (along  the  position  vector). 
East,  and  North 

h  angular  momentum 

i  inclination 

J2  second  zonal  harmonic  coefficient 

K.  E.  kinetic  energy 

k  Gaussian  or  planetary  constant 

h  additive  secular  acceleration  vector  that  accounts  for 

consideration  of  non-two-body  reference  orbit 
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SP 


true  longitude 

rotating  vectors  always  in  nominal  orbital  plane  that  maintain 
constant  angles  with  precessing  line  of  apses  of  nominal  ellipse 

mean  anomaly 

nominal  anomalistic  mean  motion 

mean  motion 

partial  double  precis -on 

potential  energy 

Kepler ian  Period 

derivative  of  Legendre  function  with  respect  to  sin  <p 

semi-parameter  of  the  conic 
rotation  matrix 
orbit-plane  system 
geocentric  distance  of  vehicle 
position  vector 

nominal  trajectory  position  vector  at  time  t 
velocity  vector 

vehicle  velocity  vector  relative  to  a  rotating  atmosphere 
acceleration  vector 

noncentral  gravitational  acceleration  vector 
atmospheric  drag  acceleration  vector 


root-sum-square 
single  precision 
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S(R,  T,  C,  R',  T',  )  state  vector  in  orbit-plane  system  coordinates 

S(r,  r,  r)  vehicle  state  vector  of  position,  velocity, 

and  acceleration 

S(r  ,  r  ,r  )  a  non-Keplerian  state  vector 

— n  — n  — n 


S(r, r,t) 
TAP 


T.E. 
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vehicle  state  vector  of  position  and  velocity 
Trajectory  Analysis  Program 
total  energy  (K.  E.  +  P,E.) 

earth-fixed  local  horizontal  system  transformation 
matrix 

time 

argument  of  latitude 
true  anomaly 
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rate  parameter  of  nominal  anomalistic  mean  motion  U 

independent  variable  difference  (time  At;  true  anomaly 
Af;  eccentric  anomaly  AE) 

Krone  eke  r  delta  function 

convergence  criterion  for  solution  of  Kepler's  equation 
vector  cross-track  direction 

rate  parameter  of  the  major  axis  of  the  nominal  ellipse 
vector  in-track  direction 
longitude  of  vehicle 

square  of  the  gravitational  constant  of  the  earth 
sum  of  celestial  object  mass  and  central  mass 
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<■  vector  radial  direction 

rate  parameter  of  the  line  of  nodes  of  the  nominal  plane 
latitude  of  the  vehicle 

longitude  of  the  ascending  node 
argument  of  perigee  or  pericenter 


SECTION  I 


INTRODUCTION 

A.  PURPOSE,  SCOPE,  AND  LIMITATIONS 

This  document  is  intended  to  serve  as  a  technical  reference  manual 
and  basic  usa^e  guide  for  the  Trajectory  Analysis  Program  (TAP).  Trajec¬ 
tory  generation  and  comparison,  as  defined  in  Section  II,  are  the  result  of 
the  evaluation  of  a  particular  formulation  for  the  equations  of  motion,  which 
yields  the  position  and  velocity  components  of  a  vehicle  at  some  desired  time 
t,  and  comparison  of  that  solution  with  one  obtained  by  another  method.  The 
mathematical  formulations  implemented  in  TAP  consist  both  of  numerical 
integration  methods  (Cowell  and  Encke)  and  analytic  algorithms  (Pines, 

Kyner,  and  Escobal). 

A  precision  numerical  integration  scheme  used  in  the  Cowell  and  Encke 
formulalions  is  described  in  Appendix  A.  A  comprehensive  description  of  all 
the  methods  used  in  TAP  is  presented  with  emphasis  placed  on  key 
characteristics. 

Sections  III  (Reference  Systems  and  Equations),  TV  (Cowell  and  Encke 
Implementation),  V  (Program  Structure  and  Logic  Flow),  and  VI  (Program 
Usage)  are  sufficiently  detailed  to  meet  the  needs  of  users  concerned  with 
trajectory  generation  and  comparison. 

Only  the  class  of  elliptical  geodetic  orbits  within  a  specific  inertial 
reference  system  is  considered;  force  models  allow  for  the  noncentral 
gravitational  effects  and  a  basic  atmospheric  drag  effect.  Particular  empha¬ 
sis  is  placed  on  the  use  of  the  Modified  Encke  formulation  for  trajectory 
computation.  The  algorithm  fo  *  the  Kyner  II  nominal  trajectory  is  developed 
in  Appendix  B. 

B.  HISTORICAL  BACKGROUND 

In  recent  years*  trajectory  generation  has  been  performed  primarily 
by  numerical  integration  of  the  equations  of  motion  in  the  Cowell  formulation 
by  the  TRACE  Orbit  Determination  Program  (Ref.  1)  at  Aerospace  Corporation. 


TAP  was  developed  to  analyze  different  computational  methods  for  trajectory 
generation,  and  to  provide  a  tool  to  establish  solution  accuracy  by  means  of 
comparison.  The  program  is  written  entirely  in  CDC  6000  series  FORTRAN 
and  its  flexible  design  permits  the  addition  of  new  formulations  and  accuracy 
tests  of  interest. 

C.  SUMMARY  OF  KEY  FEATURES 

By  a  collection  of  computational  algorithms,  TAP  generates  trajectories 
and  associated  quantities,  which  are  then  compared  and  illustrated.  The  key 
features  of  TAP  can  be  identified  by  the  following: 

1.  Methods  of  Trajectory  Generation  considered  are  the  Cowell 
formulation  for  the  reference  solution;  and  the  Encke,  Pines, 
Kyner,  or  Escobal  formulation. for  the  comparison  solution, 

2.  Event-  Detection  During  Generation  consists  of: 

•  Ascending  and  descending  nodes 

•  Apsides 

•  Specified  print  times 

•  Crash  radius 

•  Integration  step -size  change 

3.  Comparison  of  Orbital  Parameters,  both  directly  and  indirectly 
(in  the  form  of  differences  and  root- sum- square),  are  in  three 
.basic  coordinate  systems:  Earth-Centered  Inertial  (ECI), 
Classical  Elliptic,  and  Orbit-Plane. 

4.  Printer  Plot  of  prespecified  parameters  versus  time. 

5.  Integration tClo sure  results. 
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SECTION  II 


TRAJECTORY  GENERATION  AND  COMPARISON 

A.  THE  BASIC  PROBLEM 

Fundamentally,  trajectory  generation  is  a  numerical  process  for 
obtaining  the  position  and  velocity  components  of  the  vehicle  state  vector 
S(r,  r,  t)  at  some  desired  time  t.  A  basic  problem  of  trajectory  generation 
is  to  choose  a  process  that  evaluates  a  given  formulation  for  the  state  vector 
by  means  of  a  computational  algorithm.  In  the  Trajectory  Analysis  Program, 
several  mathematical  formulations  of  the  equations  of  motion  and  associated 
computational  algorithms  are  considered;  both  numerical  integration  methods 
(Cowell  and  Encke)  and  analytic  algorithms  (Pines,  Kyner,  and  Escobal)  are 
implemented. 

Trajectory  comparison  of  different  processes  or  methods  of  trajectory 
generation  is  of  particular  interest.  Typically,  the  state  vector  and  related 
quantities  resulting  from  one  method  are  compared  at  some  common  time  t 
with  those  of  another.  Numerical  properties  are  expressed  in  exact  and  dif¬ 
ference  form.  To  exhibit  a  comparison,  TAP  generates  a  precise  numerical 
integration  method  as  a  reference  trajectory,  which  is  then  differenced  at 
times  of  interest  with  Q  trajectory  generated  by  another  method.  In  addition, 
a  quick-look  plotting  feature  graphically  presents  the  behavior  characteristics 
of  specified  quantities. 

R.  TRAJECTORY  GENERATION  TECHNIQUES 

Solution  of  the  equations  of  motion  in  TAP  can  be  expressed  in  analytical 
form  by  an  algorithm  (Pines,  Kyner,  ir  Escobal),  or  it  can  be  obtained  in 
numerical  form  by  special  perturbation  methods  (Cowell  or  Encke). 

In  the  case  of  the  Cowell  and  Encke  formulations  for  numerical  inte¬ 
gration  of  the  equations  of  motion  or  deviation,  a  predictor/ corrector  (eighth- 
order  differencing)  Runge-Kutta/ Gauss -Jackson  technique  is  used  (see 
Appendix  A),  with  time  as  the  independent  variable. 
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Analytic  methods  available  in  TAP  generate  the  Keplerian  or 
non-Keplerian  state  vector  at  some  specified  time  t  by  the  evaluation  of 
closed-form  expressions,  as  follows. 


Method 

Closed-Form  Representation 

Pines 

F  &  G  functions  (Keplerian  only) 

Kyner  I 

Classical  elements 

Kyner  II 

F  &  G  functions 

E  scobal 

Classical  elements 

In  addition,  Hermitian  interpolation  techniques  (fifth-  and  third-degree 
polynomial  schemes)  are  utilized  to  obtain  trajectory  quantities  at  special 
events  (such  as  nodes,  apsides,  and  crash  radius),  and  specified 
print  times . 

C.  TRAJECTORY  COMPARISON  TECHNIQUES 

O 

The  state  vectors  and  related  quantities  at  time  t  are  obtained  from 
different  methods  and  compared  in  the  form  of  differences  (A  =  reference 
minus  comparison).  Differences  are  given  in  three  coordinate  systems  — 
Earth- Centered  Inertial  (ECI),  Classical  Elliptic,  and  Orbit-Plane  —  within 
c  which  the  maximum  separation  between  the  two  state  vectors  at  time  t  is 
expressed  as  the  root- sum- square  (rss). 

The  Cowell  formulation  has  been  chosen  to  generate  the  reference  tra¬ 
jectory  in  TAP,  thereby  defining  the  integration  step  times  as  comparison 
times.  The  comparison  trajectory  can  be  specified  as  the  result  of 

*  O 

<• 

o 

•  Encke  formulation  -  Classical  (Keplerian  nominal)  or  Modified 

'  ,  (non-Keplerian  nominal) 

•  Analytic  algorithm  -  Pines,  Kyner  I,  Kyner  II,  or  Escobal  nominals. 

That  is,  both  the  reference  and  the  comparison  trajectories  are  generated 
simultaneously  for  any  trajectory  generation  case.  In  particular,  the  follow¬ 
ing  combinations  of  methods  are  available. 

<  O 

o  j 

o  -  0  c 
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Formulation 


Reference  Trajectory  Comparison  Trajectory 

Cowell  Cowell  Analytic  Algorithm 

PINES 
KYNER  1 
KYNER 2 
ESCOBAL 

Encke  Cowell  Eneke 

CLASSICAL 

MODIFIED 

Another  useful  trajectory  comparison  technique  is  contrasting  integration 
closure  errors  obtained  from  two  different  methods.  Integration  closure 
error  is  the  magnitude  of  the  difference  vector  between  the  initial  position 
vector  £ q  and  the  resulting  position  vector  £(t  =  epoch)  after  integrating 
forward  to  some  terminal  time  and  then  backward  to  epoch.  Thus 


^CLOSURE  =|£0  "  £(t  =  epoch)l  * 


SECTION  in 


REFERENCE  SYSTEMS  AND  EQUATIONS 


The  symbols  used  throughout  this  report  follow,  after  which  the  basic 
reference  systems  and  equations  used  in  TAP  are  defined  and  discussed. 

A.  SYMBOLS  (NOTATION) 

.  time  differentation  and  vector  dot  product  operation 

. .  time  differentiation 

a  vector;  for  example,  r.  is  the  vehicle  position  vector 

(  )  relates  to  a  functional  relationship,  such  as  S(£,  r)  the  vehicle 
state  vector,  which  is  a  function  of  :r  and  r,  or  S(a,  e,  i,  $2  ,  w , 
M),the  state  vector  expressed  in  terms  of  elliptical  elements. 
Also,  used  to  denote  row  vectors;  for  example,  r  =  (x,  y,  z) 

A  increment  or  difference 

x  vector  cross  product  operation 

j  |  absolute  value,  or  vector  magnitude  operation;  for  example, 

r  =  |rj 

V  vector  gradient  operator 
||  a  matrix,  or  a  set  of  elements 
l[  ]  the  integer  operator  of  a  scalar 
[  ]  a  matrix;  for  example. 


r_  =  [x,  y,  zj  (row  vector),  or  £  =  I  y  I  (column  vector) 


1=  0 


0  3X3  identity  matrix 


L°  0  *J 


A  =  [A]  6  X  6  matrix  or  array  of  trajectory  quantities. 
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c 


computed  or  simulated  value 
i  particular  index  relative  to  a  set  or  vector 
T  matrix  transpose  symbol 

*  particular  value  of  a  predetermined  constant 


Subscripts 


C  comparison  solution  or  value 

e  the  earth 


g  the  Greenwich  meridian 

particular  index  denoting  members  of  a  set  or  components  of  a 
vector.  In  general,  i,  j,  k,  i  =  1,  2,  .  .  . 

n  nominal  value 


a  a  transformation  that  is  a  function  of  a 

y  differentiation  of  a  vector  with  respect  to  the  parameter  y 

< 

$  a  transformation  that  is  *a  function  of  $ 
o  initial  value  or  reference  point 


B.  •'  REFERENCE  SYSTEMS 

Trajectory  generation  ant.1  comparison  involve  coordinate  transformations 
and  require  selection  of  the  proper  reference  systems  for  evaluation  and  anal¬ 
ysis.  Several  such  systems  are  used  in  TAP.  Basic  to  all  is  the  earth- 
centered  inertial  {ECI)  system,  which  has  as  its  fundamental  plane  the  true 
equator  at  epoch;  its  principal  axis  is  the  mean  equinox  at  midnignt  of  the 
date  of  epoch.  This  reference  system  does  not  take  into  account  precession 
and  nutation  effects. 

1-  EARTH-CENTERED.  INERTIAL  SYSTEM  (ECI) 

Inertial  position,  velocity,  and  acceleration  vectors  are  oriented  ir.  the 
reference  system  illustrated  in  Figure  1. 

The  position  and  velocity  components  of  a  vehicle's  state  vector  within 
this  reference  system  may  also  be  expressed  in  terms  of  the  classical  elliptic 
and  orbit -plane  systems. 


radius  vector  (x,  ye  z) 
velocity  vector  (x,  y,  z) 

direction  of  the  vernal  equinox  T  in  the  equator  plane 

position  component  in  the  X  direction 

velocity  component  in  the  X  direction  (relative  to  O) 

raght-hand  axis  to  X  and  Z 

position  component  in  the  Y  direction 

velocity  component  in  the  Y  direction  (relative  to  O) 

direction  north  perpendicular  to  the  equator  plane 

position  component  in  the  Z  direction 

velocity  component  in  the  Z  direction  (relative  to  O) 

origin,  which  coincides  with  the  center  of  the  earth 

Figure  1*  Earth-Centered  Inertial  Reference  System 


3.  ORBIT  -PLANE  SYSTEM 

The  state  vector  in  the  orbit-plane  system  (Fig.  3)  is  given  at  time  t  by 

S  =  S(R,  T,  C,  R',  T',  C') 


In  this  -instance  the  point  P,  position  of  the  vehicle  in  the  ECI  system,  repre¬ 
sents  the  origin  of  the  orbit -plane  system  defined  by  the  vectors  £, 
which  are  referred  to  as  the  radial,  in -track, and  cross -track  direct^pns, 
respectively.  In  the  smaller  sketch,  it  should  be  noted  that  the  Ja  axis  is  an 
extension  of  the  geocentric  radius  vector.  The  T  axis  is  both  normal,  to  the 
j;  axis  and  positive  in  the  same  general  direction  as  the  inertial  velocity 
vector  r.  Both  lie  in  the  instantaneous  orbit  plane.  The  axis  is  normal  to 
ihe  orbit  plane  (defined  by  and  ]])  and  thus  forms  a  right-hand  orthogonal 
system. 
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The  position  and  velocity  of  an  alternate  point  Q  relative  to  point  P  are 
then  given  by 

<y 

Q  =  Q(R,T,C,R',t',C') 


where 


o  i 

R  =  position  component  in  the  j;  direction 
T  =  position  component  in  the  T  direction 

// 

\\  C  =  position  component  in  the  ^direction 
vR7  =  velocity  component  in  the  ^direction 
T7  =  velocity  component  in  the  ^direction 
C'=  velocity  component  in  the  direction 


c.:  transformation  equations 

The  transformation  of  the  state  vector  to  the  other  systems  involves 
the  following  equations. 

1.  FROM  EC1  TO  THE  CLASSICAL  ELLIPTIC  SYSTEM 

To  obtain  the  classical  elements  from  ECI  coordinates,  the  following 

O  ' 

equations  are  solved. 


S  =  r/aE  =  r  r/(|ia)^^ 

v 


<  0 


E  =  (1/ r)(PY  a) 


e2=  S2  +  C2 
e  e 


U  =  {  rh\ 

W  =  j(£  x  r>/ 1  x  r|  J 
V  =  W  x  U 


=  +  (u2  4-  V2) 

'  z  z/ 


cos  i  =  +  [(Ux  +  vy)“  +  (Uy  -  Vx) 

sin  u  =  U  /  sin  i 
z 

cos  u  =  V  /  sin  i 
z 

sin  1  =  (Uy  -  Vx)/(1  +  cos  i) 

cos  i  =  (U  -  V  )/(i  +  cos  i) 
x  y 

=  i  -  u 

C  =  p/r  -  1 
v  r 


Sv  =  r(p/p.)J 


p  =  a(l  -  e  ) 

cos  v  =  C  /e 
v 

sin  v  =  S  /e 
v 

u  =  u  -  v 


M  =  E  -  e  sin  E 
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2.  FROM  ECI  TO  THE  ORBIT-PLANE  SYSTEM 

The  orbit-plane  system  is  defined,  given  r,  £,  ECI  vectors,  by  the 
unit  vectors 

_!_  =  r/r 

£=(rxr)/|rxr| 

<1  = 


Any  ECI  vector  can  be  expressed  in  the  orbit-plane  system  (RTC)  by 
the  following  dot  products  (where  £  and  a  are  arbitrary  ECI  position  and 
velocity  vectors): 

R  =SL  *  1 

t  =  a-l 

<  c  =  a •£ 

R'=a*  I 

'  v. 

c'=i*  i 

o 

D.  TRAJECTORY  GENERATION  EQUATION^ 

Several  sets  of  trajectory  generation  equations  are  used  in  TAP  to  com¬ 
pute  the  vehicle  state  vector  S(r,r,r  }  at  some  base  time  t.  The  mathematical 
formulations  of  these  sets  are  now  described. 

o  *' 

1.  COWELL  FORMULATION 

C 

The  direct  numerical  integration  of  the  following  vector  equation  of 
motion  is  referred  to  as  the  Cowell  method: 


r  =  -y.-=^r  .  (i  =  1,2) 
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where 


r  =  [x,  y,  z] 

r  =  III 

*Li  =  Li  +l2 

=  noncentral  gravitation 

£2  —  atmospheric  drag 
2.  ENCKE  FORMULATION 

The  numerical  integration  of  the  following  deviation -vector  equation, 
which  relates  the  departure  of  the  actual  from  some  nominal  position,  is 
referred  to  as  the  Encke  method: 


P  =  r  -  r 
—  —  — n 


=  -  fi  [£/r3  -  rjrl  ]  -  L  +  fjTj  r.  (i  =  1,2) 


where 


P  =  r  -  r 
—  —  — n 


r  =fx  ,  y  ,  z  I 
— n  i  n’  J  n’  nJ 


L  =  additive  secular  acceleration  vector  that 
accounts  for  the  consideration  of  a  non- 
two -body  reference  orbit. 


The  nominal  trajectory  £n  is  the  solution  of 


/  3 

r  =-fi?/r  +L 
— n  n  n  — 

If  L  =  (),  then  the  classical  Encke  solution  results. 
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The  |i-term  in  the  above  equation  can  now  be  manipulated  into  the 
O.  K.  Smith  computational  form  (Ref.  2) 

[1+rT6][<i.n+ £>••£] -£ 

n 

where  6  -  r^/r.  Note  that  only  time  is  expressed  as  the  independent  variable 
in  this  deviation-vectoi  and  that  no  prespecified  rectification  criterion  or  pro¬ 
cedure  is  considered. 

«  3.  ANALYTIC  ALGORITHMS 


T 

r 

n. 


Analytic  algorithms  give  rise  to  closed  form  solutions  for  the  nominal 
equation  of  motion  =  -pr^/r^  +  L.  The  three  analytic  methods  of  Pines, 
Kyner,  and  Escobal  are  of  interest  here, 
a.  Pines  Method 

■O  *£?  T  r  in  '  mmm 

_  Closed -form  expressions  for  the  F  &  G  functions  are  utilized  in  the 
Pines  method  (Ref.  3)  to  evaluate  the  Keplerian  state  vector  S(.r  £n>£.n)  at 
some  specified  time  t,  given  the  initial  state  vector  SQ(r  f  i_Q)  at  initial  time 
0  (epoch)  tQ*  The  method  can  be  summarized  by  the  set  of  equations 


tn  *  £o  +  g  <‘.r0.i0) 

r  —  f(t,  r  ,  r  )  r  +  g(t.  r  ,  r  )  r 
— n  '  9  —o’  —o'  — *o  i0 

r  =  -  pr  /r3 
— d,  rr— n  n 


with  scalar  functions 


f  =  (a/ro)(cos0-l)  +  1 
g  =  t  -  [(6-sine)/nJ  =  [sin0-(0-M)]/n 

*  ,  C  ^  : 
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f  =  -  (a  n/  rrQ)  sin© 
g  =  (a/r)(cos0-l)  +  1 

and  Kepler's  equation  relating  the  time  t  to  the  angular  variable  given  by 

M  =  nt  =  0  -  sin 0(1  -  r  /a)  -  (cos 0-1)  r  •  r  /a^n 
'  o  — o  — o 

where 


a"1  =  (2/r0)  -  (V2/n) 
n^a-3'2 


r  =  (r  .  r  ) 
o  — o  — o 


1/2 


V  =  (r  *  r  ) 

o  ' — o  —o' 


1/2 


b.  Kyner  Method 

The  Kyner  method  (Ref.  4)  evaluates  a  non-Keplerian  state  vector 

S(r  ,  r  ,  r  ),  which  includes  the  first-order  secular  effects  of  the  earth's 

oblateness.  This  technique  utilizes  the  three  rate  parameters  (q,  t  ,  y  ),  in 

the  formulation  of  nominal  r  »  where 

— n 

q  is  the  rotation  rate  of  the  major  axis  of  the  nominal  ellipse 

r  is  the  rotation  rate  of  the  line  of  nodes  of  the  nominal  plane 

Y  is  defined  in  terms  of  the  nominal  anomalistic  mean  motion 

n  =  n  (1-Y),  where  n  is  the  mean  motion  at  epoch, 
o  o 
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(1)  Kyner  Formulation  I  (KYNER  1) 

The  first  forinulation  (Ref.  4)  of  the  nominal  is  given  by  the  vector 
equation 


r 


n 


[R(n,i,u)J 


that  satisfies  the  differential  equation  r_n 


2 

ur  /r  t  L  where 
r— n  n  — 


ft  =  T(f-f  )  +  ft 

'  o  o 

i  =  io 

U  =  f  +  CJ 
W  =Tl(f-f  )  +  U 

o  o 

r  =  a  (T-e^.)/(l+e  cosf) 
n  o  o  '  o  ' 

where 


(longitude  of  the  ascending  node) 

(inclination) 

(argument  of  latitude) 

(argument  of  perigee) 

(length  of  £q) 


a  =  aQ  (semi-major  axis) 

e  =  eQ  (eccentricity) 


and  the  nominal  true  anomaly  f  is  given  as  a  function  of  time  by  the  modified 
Kepler  equation 

nQ(l-Y)t  +  Mq  =  E  -  eQ  sinE  . 

O 
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Note  that 


tan{f/2)  =  [(1  +  eQ)/(l  -  eQ)J1/2  tan  (E/2) 

Mq  =  initial  value  of  the  mean  anomaly 

2  3  » 

n  a  =  H- 
o  o 

rotation  matrix  [Rj  =  [D(f2)][C(i)]  (B{u)] 

Thus  the  nominal  trajectory  r^(t)  depends  on  nine  parameters,  which  axe 

The  six  values  a,e,i,  S2,w,M 
o  o  o’  o  o  o 

The  three  arbitrary  rate  parameters  t1,t  ,y 

If  the  rate  parameters  are  set  equal  to  zero  (that  is,  if  L  =  0),  then  the 
resultant  nominal  r_  is  Keplerian. 

(2)  Kyner  Formulation  H  (KYNER  2) 

The  second  formulation  [Ref.  5]  of  the  nominal  r.  n  is  given  by  the 
vector  equation 


r  =  F  M  +  C3  N 


where 


F,  G  =  Keplerian  scalar  functions 

M,N  "  Rotating  vectors  that  are  always  in  nominal  orbital  plane 
and  maintain  constant  angles  with  the  precessing  line  of 
apses  of  the  nominal  ellipse. 


In  detail,  these  scalar  functions  are  described  as 


F  =  l/rQ(r  cos<p)  -  (r/h)o  rn  sirup 
G  =  (r/h)Q  r^  sin$> 


where 


V  =  f  ”  fo 


h 

o 


r 

— o 


Xi, 


and  the  rotating  vectors  as 

M  =  r  +  [A]  r 
—  — o  J  — o 

N  =  r  +  [A]  r 

V  —  — o  L  J  — o 

The  matrix  [Aj  can  be  represented  as  a  sum  of  constant  matrices  [A.] 
(evaluated  at  epoch)  whose  scalar  coefficients  are  slowly  varying  functions  of 
<p  and  the  rate  parameters  rj ,  t,  y.  That  is 

[A]  =  [A(?>,  t\7  t,  y,  no,  iQ»  uQ)] 

b 

[A]  a.(„)[A.] 

j=l 
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In  summation,  the  state  vector  S  (x  ,r  ,  x  )  is  expxessed  as 

v  ““H  “"H  — -n 

in  =  4+(Al]-0+G[I+[Aj]io 

£n  =-  [f[I+[A]]  +  F[A]Jxo  4  [g!I4[A]]  4  g[a]]x 

x  =  -  |iX  7x^  +  li 
— n  n  n  — 

where 


h  =  (}  -  (1-  rfWlJxl  +  (F;[A])ro  4  <Gj[A]) ^ 


(Piq) 


(differential  operator) 


where  p  and  q  a x*  the  arguments  of  the  operator, 
c.  Escobal  Method 

As  in  the  Kyner  method,  the  Escobal  method  (Kef.  6)  evaluates  a  non- 
Keplerian  state  vector  S  (x  ir  ,  j?^),  which  includes  the  first-order  secular 
\jffecta  of  the  earth’s  ohlateness.  This  technique  utilizes  the  secular  rates 
of  change  in  the  three  elements 


n  =  -f-f  £  cosic 


* _ ,3  ^3  ,,  5  .  x  n 

“=+2--2T(2-?8inio*k 


n  =  no(l  -  Y) 
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where 


i 


0 

< 


I  • 


o 

o 

<> 


Y  =  - H  (^)2(*t  (!  •  3  0iAo  S““o» 

^2  =  second  harmonic  coefficient 
Po  =  semi-parameter  of  the  orbit  a,,  epoch 
i  »  inclination  of  the  orbit  epoch 
k  =  Gaussian  or  planetary  constant 
nQ  =  k[p)  '  aQ“  '  (mean  motion  at  epoch) 

p  =  sum  of  celestial -object  mass  and  central  mass 

a  =  semi-diameter  of  the  central  mass 
e 

aQ  =  semi-major  axis  of  the  orbit  at  epoch 

rQ  =  radius  of  object  at  epoch 

uq  =  f Q  +  o  (argument  of  latitude  at  epoch) 

=  argument  of  pericenter  at  epoch 
f  =  true  anomaly  at  epoch 


(1) 


Note  that  the  formula  for  n  ie  not  the  one  given  in  Reference  6.  To 
insure  thac  the  rate  parameter  Y  for  mean  motion  is  consistent  with  the  first- 
order  secular  effects  for  SI  and  w  at  perigee,  the  secular  effect  for  the  anom¬ 
alistic  mean  motion  should  be  used;  that  is,  the  approximation  of  the 
Cunningham  integral  given  in  Eq,  (1)  above. 
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T' 


The  state  vector  S(r  n»£r,  ir  can  be  expressed  as 


=  k2  Ko£  +  2io,i  +  iiuP+yuS+2yu  Q+ywQJ 


where  the  unit  vectors 

P,  =  jP(w»  S2» i) 
Q  -  Qfw,  S2>  i) 
w  =  W(w,n,i) 

with 


p  " 

cos  w 

sin  oj 

0  ‘ 

'  1 

0 

0 

Q 

= 

-sin  w 

cos  w 

0 

0 

cos  i 

sin  i 

W 

0 

0 

i 

6 

-sin  i 

cos  i 

l  J 

' 

m 

>* 
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E.  MISCELLANEOUS  EQUATIONS 


Following 


where 


are  equations  for  specific  trajectory  quantities. 
Angular  Momentum,  h 

h  =  \lx  i| 

Root-Sum -Square,  rss 


rss  =  (R2  +  T2  +C2)1/2 


R  =  Ar.  •  £ 
T  =  Ar  •  q 
C  =  A£  •  £ 
Ar  =  r  -  r 


Keplerian  Period, 


,  3/2 

TS  -  2lra 

*k  "  .,1/2 


Kinetic  Energy,  K„  E. 


K.E.  =Yc/2 

=  r  •  r/2 
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Potential  Energy,  P.  E.  ( J ^  Only) 

P.E.  =  -  (ji/r)  [l  +  (a2  J2/2r2)(l  -  3  sin2  6)] 

where 

.2,  2,  2 

sin  o  =  2.  /r 

Total  Energy,  T.  Et 

T.  E,  =  K.  E.  +  P.  E. 

/<  1 

y  „  Gravitational  Acceleration  Due  to  the  Noncentral  Force  Field 
of  the  Earth,  r^  { J 2  Effect  Only) 


where 


U  =  “T 


i  i  .  3Z  x  * 
(_1  +-g.)-r 

r 


a  =  1  earth  radius 


Noncentral  Earth  Gravitational  Force,  r  , 

_ _ 

The  gravitational  acceleration  due  to  the  noncentral  force  field  of 
the  earth  is  derived  from  the  generalized  potential  function 

n 

U  =  (p/rf£  (%/r)n^  P™  (sin  ?>(cnm  cos  m  \  +  Snm  sin  m  x) 
n=2  m=o 

where 

r,  cp,  X  are  the  geocentric  distance,  latitude,  and  East 
longitude  of  the  vehicle 

ag  is  the  mean  equatorial  radius  of  the  earth 

Pm  is  the  Legendre  associated  function  of  the  first  kind 
n  of  degree  n  and  order  m 

C  ,  S  are  numerical  coefficients, 
nm  nm 

Components  of  the  acceleration  vector  r^  are  expressed  as 
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■v 


where 


Note  that 


so  that 


T  is  the  earth-fixed  local  horizontal  system 
<P’  transformation  matrix 

— LH  *s  t^ie  gravitati°nal  acceleration  due  to  the  noncentral 
force  field  of  the  earth  in  the  local  horizontal  system, 
in  which  the  coordinate  axes  are  directed  Up  (along 
the  position  vector),  East,  and  North. 


-EH  =  °East 


North- 


w 

Gup  =3U/3r  K*  (f.M 

n=2 

CO 

GEast  =  (1/r  C°S  Kn  (<P'X) 


where 


■w 

•  SNorth=<1/r><8U/W>=EKn(<,’X) 


(¥>, X.)  =  -  (n/r2Hn+l)(ae/r)n£]  P™(ain?>)  toimH  sin  m  \] 
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>  * 

> , 


n 


K2  (<p,  X)  =  (n/r2  cos  v>)(ae/r)n^^  m  P“  (sin  <p)  jc 


sin  m  X  -  S  cos  m  x] 
nm  nm  J 


tn=o 


Xu 


K2  (?,X)  =  ^/r2)(ae/r)n^  P™  (sin  <p)  cos  pjc^  cos  m  X  +  Snm  sin  m  x] 


m=o 


(Pm  is  the  derivative  of  the  Legendre  function  with  respect  to  sin 9.) 


n 


Two  different  optional  normalizations  of  the  Legendre  functions  are 


used.  The  alternate  values  for  C  and  S  coefficients  (Ref,  7)  are  denoted 

nm  nm 

by 


C  ,  S  (APL  normalization) 
nm’  nm 


and 


where 


and 


C  ,  S  (Kaula  normalization) 
nm  nm 


(Cnm-  5nm)  =  [(n+n,)!/(n-m)!]1/2  (C^,  SnJ 


i 

% 

C, 


l  ' 


(C_t  S.„)  =  |" (En-l)(2-6pm)]"1/2  F_w_j 


AAAAA  AAA117  U 


-*  *  imi~  mu* 


with  6.j  being  the  Kronecker  delta  function. 
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a'if  ■'ill  I 


A 


Atmospheric  Drag  Acceleration,  r  ^ 


r2  =  - P(va/2)(cda/W) 


•where 

v 

3  1  K 

p  =  the  density  at  height  h  above  an 
oblate  earth  (the  roily  atmosphere 
model  considered  is  the  Lockheed- 
Jacchia) 

O 

CDA/W  =  the  drag  coefficient 

=  vehicle  velocity  vector  relative 
to  a  rotating  atmosphere 

>  '  ~A  =  [XA’  ^A*  2a] 


where, 


c* 


o 


v 


*A  =  *  +  “la  y 
yA  *  y  -  “!a  x 


<*>a  =  rotation  rate  of  the  atmosphere. 
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Apsis  Event  Detection  (between  times  t.  and  t._j) 


If  cosp^  ♦  cosp.  ^  >  0,  there  is  no  event.  Otherwise,  an  event  has 
occurred  and 


9  cosp./ 8 1  =  r .  .  r  .  +  r  .  •  r  .  for  all  i  . 
ri  — 1  —x  —x  —1 


Since  cosp  is  a  smooth  monotonic  function  of  t  and  cos  p  is  near  zero,  then 


at/3  cosPj^  =  (3  cospj/at) 


-1 


and  allows  interpolation  for  the  time  t  for  which  cos  p  =  0,  given 


cosp.,  cosp^,  t.,  U  v 
8t./8  cospp  8t.  j/3  cosp..  ^ 

After  tCQ8  p_Q  is  obtained,  interpolation  for  the  state  vector  yields 

S(r,r;  tcos|J=0> 

given  tuv  r.,  r.,  r.,,  r., 

Node  Event  Detection 

If  z^  .  2  0,  there  is  no  event;  otherwise,  an  event  has  occurred. 

Then  assume  that 


8ti/8zi  =  (Sz^/etj) 


-1 


-31- 


i 


and  interpolate  for  the  time  at  which  z  ~  0.  After  is  obtained, 
interpolation  for  the  state  vector  yields 


i;  W 


SECTION  IV 


COWELL  AND  ENCKE  IMPLEMENTATION 


Table  I  indicates  the  sequence  of  steps  necessary  to  implement  Cowell 
and  Encke  formulations. 


Table  I.  Schematic  for  Implementing  Cowell  and  Encke  Formulation* 


Note  1:  Subscript  n  denotes  nominal  values  obtained  by  some  analytic  method. 

Note  E:  If  the  user  chooses  to  specify  rectification  criteria  e  and  i,  then 
investigate  the  rectification  criterion  IJPja  «,  or  |£| a  «.  If  rectifi¬ 
cation  necessary,  set  SQ(r  q,  r  .)  =  S^r.r)  and  proceed  with  Step  2. 
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SECTION,  V 


PROGRAM  STRUCTURE  AND  LOGIC  FLOW 


The  general  structure  and  logic  flow  are  illustrated  in  Figure  4. 


*TAP* 


I 


«.  ■ 
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MULTIPLE 

CASE 

LOOP 


(STORAGE) 

(DATA) 

[preset  region] 
-“[general  input  region] 

•  input  data 

V 

[iNITIALIZ  ATIONsREGIOn] 

•  INITIALIZE  FILES 
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Figure  4.  Diagram  of  TAP  General  Structure  and  Logic  Flow 


SECTION  VI 


PROGRAM  USAGE 


The  general  use  of  the  Trajectory  Analysis  Program  will  be  shewn  by 
describing  the  specific  input  data,  output  tables,  and  plotting  features  avail¬ 
able  to  the  user. 

A.  INPUT  DATA 

All  required  input  is  nominally  set  to  values  in  the  basic  unit  system  of 
feet,  degrees,  and  feet  per  second.  In  Table  II,  each  input  variable  is 
described  and  its  preset  value  designated  in  the  value  field. 


Table  II.  Input  Data  Description  (cont.) 


Code 

Location 

Value 

Description 

I 

IPLX&T 

0 

Jplot  file  data  generation  flag 

IPL&T 

=  o  , 

No  plot  file 

-  n  , 

Plot  file  generated 
every  n  steps 

I 

IUNIT 

1 

Initial  units  conversion  flag 

<*> 

IUNIT 

=  o  , 

No  conversion 

'' 

=  1  , 

Convert  from  ex¬ 
ternal  units  to  in¬ 
ternal  units 

I 

IPRNT 

r 

Print  flag 

IPRNT 

=  0  , 

Only  difference  print 

=  1  , 

Standard  print 

U 

PSTEP 

60. 

Print  step 

PSTEP 

*  0  , 

Trajectory  print 
every  PSTEP  inter¬ 
val  from  TZERO 

I 

ID  IFF 

0 

Integration  closure  flag 

IDIFF 

=  0  . 

Closure  not 
considered 

£  0  . 

Closure  error 
computed 

TZER<& 

(*> 

Reference  time  (epoch) 

TEND 

(*) 

Termination  time  relative  to  TZER0  =  0. 

:  XO 

<*> 

X 

o 

2 

(*) 

yQ  ! 

Initial  position  vector 

t 

3 

(*) 

7. 

0 

_ 

( 

v*r\n 

2 

111 

i 

> 

yo 

i 

Initial  velocity  vector 

‘ 

3 

I 

z 

o 

c 

- 

'"Required  for  each  trajectory  case. 


Table  II.  Input  Dwta  Description  (cont.): 


Code 

Location 

Value 

Description 

•  •  • 

CRASH 

1. 

Crash  radius 

GM 

.55303935E-2 

Gravitation  constant 

CJ2 

1082. 3E-6 

J2 

EPS 

1.0E-9 

Analytic  convergence  criterion 

DRAD 

(57. 295779513082) 

Degrees  per  radian  constant 

CPI 

(3.14159265358979) 

ir 

♦  •  • 

•  •  • 

DF 

20925738. 

Distance  conversion  factor 

> 

VF 

248762.3 

Velocity  conversion  factor 

AF 

5812.705 

Acceleration  conversion  factor 

TF 

60 . 

Time  conversion  factor 

HO 

1. 

Initial  integration  step 

HMIN 

.015625 

Minimum  allowed  integration 
step 

HMAX 

64 . 

Maximum  allowed  integration 
step 

ER 

L0E-10 

Integration  truncation  control 

I 

IR 

8 

Ratio  of  Runge-Kutta  steps  to 
Cowell  steps 

TDJO 

0. 

Julian  date  of  epoch. 

D1 

6.83  \ 

Coefficients  used  in  the  Lockr 

D2 

-15.684  j 

heed-Jacchia  atmosphere  model 
for  the  density  expression. 

FLUX 

0. 

10.7  cm  solar  radiation  in  units 
of  10"20  watt/m^;  if  equal  to 
zero,  the  FLUX  is  computed 
internally. 
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ft.  - 


II g—- 


Code 


Table  II.  Input  Data  Description  (cont.) 


Location 


Value 

Description 

(*) 

C][>A/W:  Cjj  is  aerodynamic 
drag  coefficient;  A  is  the  aver¬ 
age  cross-section  area  of 
vehicle;  W  is  the  weight  of 
vehicle;  units  are  ft^/lb. 

If  option  used  input  required  for  each  trajectory  case 


OUTPUT  DATA 


Tables  III  through  IX  describe  the  basic  output  quantities. 


Table  III.  Initial  Condition  Data  (Internal  Units) 


X 

x  1 

F„ 

TZERO 

o 

° 

xo 

o 

o 

M 

-  4  FVo 

F*  TEND 

—  o 

z 

* 

Fz 

GM 

o 

O 

zo 

J 

J 

DF 

VF 

AF 

TF 

❖ 

Not  available 


Table  IV.  Epoch  (TZERO)  Data  (Internal  Units) 


T 

X 

* 

X 

NSTEP 

y 

l  y 

» 

£ 

•  • 
y 

•  • 

£ 

HO 

z 

• 

z 

•  • 
z 

H  (current  integration 
step) 

XC 

*c 

••  1 
X. 

N  (#  of  equations) 

yc 

lc  yc 

•  ic 

•  • 

y 

•  • 

•  21c 

zc 

•  • 
z 

d 

h 

X 

Ax 

Ax 

Ax 

• 

h 

y 

.  h  =  £  X  £ 

Ay 

>  Ar  =  r-r c  Ay 

.  A£ 

Ay 

c 

i  •• 

’  AE 

h 

z 

Az 

Az 

,  Az 

h 

(|h|) 

d 

d 

Table  IX.  WKCNC2  Package  Output  (Internal  Units) 


C.  PLOTTING  FEATURES 


Printer  plots  are  given  for  prespecified  orbital  parameters  versus 
time.  The  parameters  considered  are: 

RSS  Root-sum-square  of  difference  between  Reference  and 
Comparison 

ENERGY  Difference  in  total  energy  {TE  =  KE  +  PE)  between 
Reference  and  Comparison 

§(XI)  Radr'l  component  of  difference  between  Reference  and 
Comparison 

q  In-track  component  of  difference  between  Reference  and 

Comparison 

r(R)  Radius  of  Reference 

a(R)  Ssmimajor  axis  of  Reference 

e{R)  Eccentricity  of  Reference 

i(R)  Inclination  of  Reference 

S2(R)  Right  ascension  of  the  ascending  mode  of  Reference 

w(R)  Argument  of  perigee  of  Reference 

M(R)  Mean  anomaly  of  Reference 

TE(R).  Total  energy  of  Reference 

hz(R)  Polar  component  of  angular  momentum  of  Reference 

r(C)  Radius  of  Comparison 

a(C)  Semimajor  axis  of  Comparison 

e{C)  Eccentricity  of  Comparison 

i(G)  Inclination  of  Comparison 

£2(C)  Right  ascension  of  the  ascending  mode  of  Comparison 

w(C)  Argument  of  perigee  of  Comparison 

M(G)  Mean  anomaly  of  Comparison 

TE{C)  Total  energy  of  Comparison 

hz{C)  Polar  component  of  angular  momentum  of  Comparison 

i.i>iuvv»u  ouci  gy  ui  uuuijjcix'iouil 

Examples  of  a  subset  of  these  printer  plots  are  given  in  Figures  5 
through  & . 
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Figure  6.  Printer  Plot  of  Variable  =  (XI) 
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Figure  8  .  Printer  Plot  of  Variable  =  M,  C 
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NUMERICAL  ACCURACY  OF  THE  DEQ  NUMERICAL  INTEGRATOR 

1.  INTRODUCTION 

In  recent  months,  much  experience  has  been  gained  at  the  Aerospace 
Mathematics  and  Computation  Center  with  a  particular  numerical  integrator 
(Ref, 8),  referred  to  as  DEQ  and  used  principally  for  trajectory  generation. 

Developed  and  written  by  J  ames  F.  Holt  of  Aerospace,  it  is  designed  to 
numerically  integrate  a  set  of  N  simultaneous  2nd-order  ordinary  differential 
equations.  A  fourth-order  Runge-Kutta  method  is  used  for  its  starting  pro¬ 
cedure  and  a  Gauss -Jackson  method  (with  eighth -order  differences)  is  use . 
for  normal  integration.  Additional  technique!  gllow  a  variable -step  mode 
(halving  and  doubling)  with  automatic  local  truncation-error  controls . 

In  trajectory  generation,  DEQ  is  used  to  integrate  the  equations  of 
motion  expressed  in  a  Cowell  and/or  Encke  formulation.  The  numerical 
accuracy  of  this  integrator  for  typical  geodetic  orbits  is  of  primary  interest 
here.  To  investigate  the  solution  accuracy,  direct  comparisons  were  made 
with  known  analytic  solutions  and  integration  closure  tests  were  performed 
with  augmented  force  models. 

2.  BASIC  MODEL 

A  Cowell  formulation  of  the  equations  of  motion  in  vector  form  was 
used  for  testing  where  the  total  acceleration  vector  acting  on  a  vehicle  was 
represented  as  the  sum  of  a  primary  gravitational  term  and  total  perturba-  »  ! 

tive  term.  Table  A -I  shows  such  a  vector  in  an  ECI  coordinate  system.  V 
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Problem: 


Numerically  integrate  second-order  ordinary  differential 
equations  of  the  form 


-  mV-L 


+  L  li 
i-l  1 


where 


N  . 

Lv.M 

i=l~1 


-  Position  vector 

=  Newtonian  Gravitational  Constant 

=  Acceleration  vector  resulting  from 
perturbing  forces 


Solution: 


DEQ  = 


Floating-point  Cowell  (second  sum)  Runge-Kutta 
integration  of  second-order  equations 

•  4th-order  Runge-Kutta  method  to  start  and 
halve  the  step- size  during  the  integration 

<*  Cowell  "second-sum"  (Gauss-Jackson)  method 
based  on  8th  differences  to  continue  the  integration. 


a.  NUMERICAL  ACCURACY 

All  numerical  results  presented  here  are  as  of  12  June  1967. 

An  example  of  the  numerical  accuracy  of  this  integrator  on  three 
different  computers  -  IBM  7094,  CPC  3600,  and  CDC  6600  -  is  given  in 
Figure  A-l.  Note  the  following  respective  word  lengths  (DP  =  double  preci¬ 


sion): 


Machine 


IBM  7094 


r* T\r*  o4nn 

v*/v  «/wv 


CDC  6600 


Word  Length 


36  bits 


ASX 


Significant  Digits 
8  (16  DP) 


i  i  toe  rvm 

x  a  f 


60  bits 


t  - 


The  root-sum-square  (rss  of  the  error)  is  the  difference  obtained  from  com¬ 
paring  a  known  analytic  solution  at  each  integration  step  and  is  plotted  vs 
time  by  using  a  log-log  scale.  Note  that  a  fixed  step  of  one  minute,  using 
6th-order  differences  with  partial  double  precision  (PDP)  on  the  IBM  7094,  gives 
a  100-ft  error  after  one  day;  whereas  the  8th-order  method  using  single  preci¬ 
sion  (SP)  on  the  CDC  6600,  gives  an  approximate  1/ 100-ft  error,  and  only  2  ft 
after  10  days.  Comparable  results  were  obtained  by  using  DP  on  the  CDC  3600. 

b.  TEST  CASES 

The  test  cases  chosen  for  this  study  are  listed  in  Table  A -II. 


Table  A-II.  Test  Orbits 


Orbit  A 

Orbit  B 

Orbit  C 

a 

21533095. 

86767518. 

21654273. 

(ft) 

.0000117 

.737 

.03115 

i 

45. 

63.4 

95.3 

(deg) 

n 

0. 

177. 

351.6 

(deg) 

U) 

H-* 

00 

o 

270. 

295. 

(deg) 

88.2 

713.4 

88.9 

(min) 

h 

P 

99.9 

317.7 

18.3 

(n  mi) 

h 

a 

100. 

21373. 

240.3 

(n  mi) 

A  parametric  relation  between  the  error  variation  rss  and  the  fixed-step  size 
h  for  a  particular  time  period  can  be  utilized  to  investigate  optimum  step 
size,  in  Figure  A-2{a)  this  technique  is  illustrated  by  plotting  the  relative 
maximum  rss  error  vs  the  integration-step  size  for  orbits  A,  B,  and  C. 

Note  the  roundoff  error  growth  for  Orbit  B  to  the  left  of  the  optimum  step 
(0.5  min)  as  h  tends  toward  zero.  On  the  other  hand,  the  high  eccentricity 
of  this  orbit  causes  sudden  changes  in  acceleration  near  perigee  and  a  very 
uniform  acceleration  at  apogee,  which  results  in  a  greater  growth  for 
h  than  0 .  5  minute . 
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Optimum  steps  for  Orbits  A  and  C  are  1.5  and  1.0  minutes,  respectively. 
Thus,  if  the  investigation  time  period  is  changed,  the  relative  optimum  step 
may  also  change;  for  example,  Orbit  C  has  a  1.5-minute  optimum  for  a  1-day 
sample . 

(>1 )  Variable  Step  Mode 

In  general,  the  variable  *  step  mode  reduces  the  rss  error  (for  a  10-day 
period).  For  example,  for  Orbit  B  in  Figure  A-2  (b),with  hQ=  0.5,  it  de¬ 
creases  the  rss;  and  for  h  =  1,  it  decreases  the  error  from  25  to  8  ft.  The 

o 

behavior  is  the  a' \me  at  the  end  of  1  day. 

(2)  Integration  Closure  Tests 

With  this  numerical  integration  method,  closure  tests  were  made  on 
three  augmented  force  models.  A  generalised  potential  function  expansion 
for  the  aspherical  earth  was  assumed,  and  the  symbols  Z,  T,  and  R  denoted 
zonal  harmonics,  zonal-plus -tesseral  harmonics,  and  harmonic s -plus - 
resonance,  respectively.  Relationships  between  closure  error  and  fixed-step 
size  permitted  the  desired  "best'  step  size  to  be  chosen  for  a  particular  force 
model  and  orbit.  Figure  A-3  indicates  that  for  a  l -day  sample  interval, 
h  =  1.0  minute  is  optimum  for  the  Z  and  T  models,  and  h  =  0.5  for  the  R 
model  at  the  desired  level  of  accuracy.  In  another  example  (Figure  A-4  for 
Orbit  C),  h  =  0.5  is  necessary  for  both  the  T  and  R  force  models. 

(3 )  Constant  Energy  Considerations 

For  simple  conservative  fields,  such  as  spherical  and  Z  models,  both 
the  magnitude  of  angular  momentum  and  its  polar  component  exhibited  linear 
growth  with  time,  but  maintained  an  acceptable  degree  of  significance  as 
defined  by  computer  word  length. 

(4)  Sample  Running  Times 

Table  A -III  gives  sample  running  times  for  trajectory  generation  by 
ufcing  three  augmented  force  models.  All  stated  times  are  considered  upper 
bounds . 
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APPENDIX  B 

KYNER  II  NON-TWO -BODY  NOMINAL  TRAJECTORY 


1.  FORMULATION 

The  formulation  of  the  Kyner  II  non-two-body  nominal  trajectory  is 
developed  in  Reference  5. 

2.  COMPUTATIONAL  ALGORITHM  FOR  THE  REFERENCE  TRAJECTORY 

To  utilize  the  theoretical  concepts  presented  in  Section  III.D.  a 
computational  algorithm  is  developed  to  evaluate  the  reference  trajectory. 

The  algorithm  is  composed  of  basic  equations  of  the  formulation,  related 
formulas,  and  the  computational  scheme. 

a.  BASIC  EQUATIONS  OF  THE  FORMULATION 

Given  the  initial  position  rQ  and  velocity  vectors  at  epoch,  the 
reference  r  can  be  written  in  terms  of  the  F  and  G  functions  as 


where 


r  =  F  M  +  G  N 


M  =£o  +  H  =  [l  +  [a]] 

*=4  +  [AK  =  MAIU 


(B-l) 


F  -  (r/rQ)  cos  $  -  (r/h)0  r  sin$ 
G  =  {r/h)Q  r  sin$> 


where 


-  v  -  v 


8 

K  («  K] 
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a^  are  slowly  varying  scalar  functions  of  $ 

[A^J  are  constant  matrices  evaluated  at  epoch 

Thus  Eq.  (B-l)  becomes 

x  =  F[l  +  [A]] ^  +  G  [i  +  lAj]^  (B-2) 

The  reference  velocity  £  is  then  obtained  by  differentiating  Eq.  (B-2)  as 
x  =  [F[I  +  LA3]+  F[A]]  x^  +  [6 [I  +  [Aj]+  GtA]]^  (B-3) 

where 

F  =  (r/r o)(ir  cos  $  -  r  sin  $$)  -  (r/n)o  (*  sin  $  +  r  cos  $$) 

G  =  ( r/h)o  (r  sin  $  +  r  cos  $$) 

r 

r  =  a  AE  |Ce  sin  A  E  +  Sg  cos  A  eJ 

AE  =  ha/  r 
8 


The  differential  equation  satisfied  by  r  then  becomes 


where 


L.*  [l  “  (1  -  V)2]pr/r3  +  (f;  £Aj)  +  (g;  [Aj)  ^ 


where 

y 


-i-m: » 


-  .  2  .  .2  . 
3  am  l  sin  u  ) 
o  o 


(p;q)  =  2  +  p  (differential  operator) 


and  p  and  q  are  arguments  of  the  operator. 

b.  RELATED  FORMULAS 

(1)  Kepler1  s  Equation  in  Differenced  Form 

The  independent  variable  is  expressed  in  difference  form  A  =  t  -  tQ 
(time),  A=  v  -  v  (true  anomaly),  or  A  =  E  -  Eq  (eccentric  anomaly).  It  is 
therefore  convenient  to  express  the  fundamental  angle-time  relationship 
(Kepler’s  equation)  as 

M-M  =(E-E  )  +  2S  sin2  (E  -  E  )/2  -  C  sin  (E  -  E  )  (B-5) 

o  o  e  o  e  o 


where 

M-M  =  n  (t  -  t  ) 
o  '  o' 

C  t  (e  cos  E) 
e  o 

S  =  (e  sin  E) 
e  o 
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( 2 )  True  Anomaly  and  Eccentric  An omaly  Relationships 

Difference  forms  relating  the  true  anomaly  and  the  eccentric  anomaly 
are  necessary  to  evaluate  r,  i_,  and  r.  These  are 


2  2 

cos  (v  -  vQ)  =  1  -  2-4Ll£  )  [l  -  COS  (E  -  Eo)] 


2  2  2 

sin  (v  -  vq)  =  -■  e  ]  [(M  -  M0)  -  (E  -  Eq)  +  sin  (E  -  Eq)] 

sin  (F<  -  Eq)  =  — sin  (v  -  vq)  -  ^  [l  -  cos  (v  -  vq)]  Se 


sin  (V  -  vo) 


cos  (E  -  E  )  =  1 - £ 

o'  ap 


1  -  cos  (v  -  vo)J 


(3)  Expressions  for  rt  r  and  <$ 

If  A  =  t  -  tQ,  or  A  =  E  -  E^,  then  the  reference  position  distance  is 
given  in  difference  form  by 

r  =  (l  -  C  cos  (E  -  E  )  +  S  sin  (E  -  E  )J  (B 

L  e  o  e  o J 


If  A  =  v  -  v  ,  then 
o 


where 


r  =  p  [l  +  Cv  cos  {v  -  vq)  -  Sv  sin  (v  -  v^)]"1 


C  =  (e  cos  v) 
v  'o 


S  -  (e  sin  v) 
v  ' o 


Differentiation  of  Eq.  (B-6)  gives 

r  =  a  A  E  Fc  sin  (E  -  E  )  t  S  cos  {E  -  E  )1 


where 


AE  =  n  (a/r) 


If  $  =  v  -  v  ,  then 
o 


$  =  h/r2 


where 

i 

h  =  n  a2  (1  ~  e2)2  =  (n/nQ) 

c.  COMPUTATIONAL  SCHEME 

If  the  initial  position  vector  2^  and  velocity  vector  are  given  at 

epoch  time  tQ,  the  following  procedure  will  establish  a  means  of  determining 

the  reference  position  r,  velocity  r,  and  acceleration  jr  vectors  at  time  t. 

# 

Two  computational  phases,  initialization  and  evaluation,  are  presented  in 
detail* 

(1)  Initialization  Phase 

All  initial  orbital  parameters,  rate  parameters,  coefficients,  and  time 
invariant  matrices  are  evaluated  at  epoch  by  the  iollowing  equations. 

(a)  Initial  Orbital  Parameters 

The  desired  set  of  orbital  parameters  is 

(r  ,  a,  C  ,  S  ,  e,  1,  u,si,C,S,  v  ^  w  ^ 
o  e  e  o  o  v  v  o  o 


<>• 


1  4 


>  4 


o 
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*  «  * 
>  « 


'v'S'F' 


where 


rQ  is  the  geocentric  radius  distance  at  epoch 

a  is  semimajor  diameter  of  the  conic  at  epoch 

C  ,  S  are  coefficients  evaluated  at  epoch 

e  is  the  eccentricity  of  the  conic  at  epoch 

i  is  the  inclination  of  the  conic  at  epoch 

uQ  is  the  argument  of  latitude  at  epoch 

£2q  is  the  longitude  of  the  ascending  node  at  epoch 

Cv,  Sv  are  coefficients  evaluated  at  epoch 

vq  is  the  true  anomaly  at  epoch 

gjo  is  the  argument  of  perigee  at  epoch 

Equations  for  the  orbital  parameters  are 
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K,  (!'«;>/ a)" 


S  r  /a  £  ~  (e  sin  E) 

O  O  O  £ 

e  =  (C2  +  S2)2 
e  e 

U  =  r  /r 
— o  — o  o 

W  =  (r  X  r  )/|r  X  r  | 

— c  — o  — o  l—o  — ol 


y  =w  x  u 

— o  — c  — o 


Bin  i  =  [u2  +  Y2  l2 
L  oz  ozJ 


COB  i  =  f(U  +  V  )2  +  ;(U  -  V  )212  -  1 

L'  ox  oy  "  o,«  ox'  J 


i  =tan  *  (sin  i/cos  i) 


sin  u 

=  U  /  sin  i 

0 

02 

COB  U 

=  V  /  sin  i 

O 

oz 

U 

rt 

=  tan  ^  (sin 

sin  i  =  (U  +  V  )A(1  +  cos  i) 
o  ox  oy 

cob  i  -  (U  -  Y  )/  (1  i  cob  i) 
o  oy  ox 

i  =  tan  ^  (sin  i  /cos  I  ) 
o  o  o 


—  j» 

“*  TO. 

o  o  o 


.1 - ..«x. .  ^  .V 

Tti  uc  xxxiu.^: 7 


p  =  a(l  -  e  ) 

C  =  p/r  -  1  =  (e  cos  v) 


(semi-parameter) 
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Sv~*Ao(p/fi)  =(esil*v)o 

siln  v  =  S  /e 
o  v 

cos  v  =  C  /e 
o  v 

v  =  tan  *  (sin  v  /cos  v  ) 
o-  o  o 

If  e  =  0,  set  v  =0.  Then 
o 


,b)  Initial  Rate  Parameters 

The  initial  value  of  the-  rate  parameters  discussed  in  Reference  4  are 
evaluated  from  the  equations 

V  =  -  |l2  (ae/a)^  (a/r)2  {1  -  3  sir2  iQ  sin2  uo) 

’i=|j2{ae/a)2(l-e2r2(4-5»in2y 

T=-|j2  (ae/a)2  (1  -  e2)'2  cos  i« 
h  =  n(l  -  Y) 

where 


-S"7”- 


c.  F  and  G  Coefficients 

The  F  and  G  functions  have  initial  coefficients  F^,  F^,  and  G^,  which 
are  evaluated  at  epoch  as 


it! 


F1  =  (1/r>o 


f2  -  -  (f /h)o 


G1  *  lr/h>o 


d.  Time  Invariant  Matrices 

The  constant  orientation  matrices  (A.)  (j  =  1,  ,  8)  discussed  in 

3 

Section  III.  D.  3.  b  are  given  as 


[Aj]  = 


0-10 


[A,]- 


0 

0 

0 

1 

0 


0 

0 

0 

0 

0 


[a3]  =  [D(n0)]  [C(io)J  [A2]  [C(-io)]  [D(-S2o)] 
[A4J  =  [D{no)J  [c(io)J  [a,]  [c(-io)]  [D(-no)j 


wnere 


[D(S2)j  = 


cosn  -sinS2  0 
sin  cosS2  0 

.  e  o  i 


i; 


O' 

$■  < 
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r  k‘ 


10  0  ' 
[C(i)J  -  0  cos  i  -sin  i 

0  sin  i  cos  i 


[A5]  =  [A  J  [A3] 
=  !A3]  ^4! 
[A7]  =  [A2]  Ia3] 


[Aq]  =  [Aj]  [A4J 


(2)  Evaluation  Phase 


To  obtain  the  nominal  state  S  (r  ,  r  ,  r  ,  L)  for  some  given  value  of 

— n  — n  — n  — • 

the  independent  variable  (time  t,  true  anomaly  v,  or  eccentric  anomaly  E), 
the  following  equations  are  solved. 

First  determine  the  value  of  the  true  anomaly  difference 


$  =  v  -  v  . 
o  o 


If  time  is  the  independent  variable  (A  =  elapsed  time  from  epoch), 


then  M  -  MQ  “  n  A ,  and  Kepler's  equation 


_  /E  -  E  \ 

M  -  M  =(E-E)  +  2S  sin  l — =r~2)  - 
o  o  e  '  2  / 


C  sin  (E  -  E  ) 
e  o 


is  solved  for  E  -  E  .  Thus 
o 


r  =  a  (l  -  cos  (E  -  En)  +  Sp  sin  (E  -  Eo)) 


(B-7) 


cos  $  =  1 


2 II  2\ 

a  (1  -  e  ) 


-  cos  (E  -  Eq)  j 


i 


(<M  -  Mo)  -  (E  -  Eq)  +  sin  (E  -  E)) 


sin  $ 


»-e¥/2 


* 


tan 


/  sin$  \ 
\cos$/‘ 


Now  apply  a  moncfonicity  adjustment  on  S>  (see  Section  B.  3.f  of  this 

Appendix).  If  true  anomaly  is  the  independent  variable  (A  »  $  &  v  _  v  )  then 
.  o 

compute  eccentric  anomaly  difference  E  -  Eq  from 


sin  (E-Eq)  =  r/(ap)1/2 


sin  $  -  (r/p)  (1  -  cos  $ )  Se 


cos  (E  -  Eq)  =  1  -  rrQ  (ap)" 1  (1  -  cos  $  ) 


where 


r  =  p/(  1  +  Cv  cos  $  -  Sv  sin  $ ) 


Thus 


M  -  Mq  =  (E  -  Eq)  +  2Sg  sin2  [(E  -  Eq)/2  J  -  Ce  cos  (E  -  E  ) 


Now  apply  a  monotonicity  adjustment  on  M  -  Mq  (see  Section  B.3.f)  and 
compute  the  corresponding  elapsed  time  from  epoch  as 


A  —  —  iXJt  x  r  \  Jzr 

*—  v  ~  \xvx  “  ™  n 11  • 
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If  eccentric  anomaly  is  the  independent  variable  (A  =  E  -  E^),  then  compute 

M  -  M  and  t  -  t  as  above  and  proceed  with  Eq.  (B-7). 
o  o 


Evaluate  the  F  &  G  functions  and  scalar  coefficients  a.  of  the  rotation 

3 


matrix  [A]  as  follows. 


F  =  r/r  cos  $  -  (r/n)  r  sin  $ 
o  o 


G  =  (r/h)Q  *  sin  $ 


F  =  (r/r  ) (r  cos  $  -  r  sin  <$$)  -  (r/h)Q  (r  sin  $  +  r  cos  $6) 


6  =  (r/h)Q  (r  sin  $  +  r  co&  $$) 


where 


r  =  a  (Alb )  sin  (E  -  Eq)  +  Sg  cos  (E  -  Eq)) 


(AE)  =  na/r 


$  =  h/r  h  =  (1  -y)h 


8  •  8 

and  {a.}  and  {a.}  ,  are  as  given  in  Eq.  (B-8)  in  Section  B.3.f. 

j  J  J  J  * 


iA]  =  £  a.  [A.] 

j  =  l 
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and 
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SUBROUTINE  WKCNC2 


a.  IDENTIFICATION 

WKCNC2:  Computational  Algorithm  for  Generating  Nominal  Trajectory 
6600  FORTRAN 
Aerospace  Corporation 

b.  PURPOSE 

To  generate  a  nominal  trajectory  {  r,  £,  £}  at  time  t  using  a  modified 


F  and  G  formulation. 

c.  NOTATION 

r 

—  o 

(XO) 

initial  geocentric  position  vector  of  the 
vehicle  at  epoch 

• 

—  o 

(XDO) 

initial  geocentric  velocity  vector  of  the 
vehicle  at  epoch 

H- 

(GRAY) 

square  of  the  gravitational  constant  of  the 
earth 

J2 

(CJ2) 

second  zonal  harmonic  coefficient 

Tt 

(PI) 

numerical  constant  (v  =  3.  1415927. . .) 

Ck 

(EFSK) 

convergence  criterion  for  solution  of 
Kepler’s  Equation 

A 

(ARG) 

independent  variable  difference  (time  At, 
true  anomaly  Af,  or  eccentric  anomaly 
AE) 

{*.  £,  r  } 

r, 

(XC) 

nominal  trajectory  point  on  the  conic  at 
time  t  where 

K 

_r  is  the  position  vector 

r  is  the  velocity  vector 

r_  is  the  acceleration  vector 

2t r 

(TW0PI) 

numerical  constant 

ro 

(RO) 

geocentric  radius  at  epoch 

1 


a 

(CA) 

semimajcr  axis  of  the  conic  at  epoch 

e 

(CE) 

eccentricity  of  the  conic  at  epoch 

i 

(Cl) 

inclination  at  epoch 

u 

o 

(CUO) 

argument  of  latitude  at  epoch 

^o 

(C20) 

true  longitude  at  epoch 

(C0) 

longitude  of  the  ascending  node  at  epoch 

V 

o 

(VOARG) 

true  anomaly  at  epoch 

0) 

(ARGPER) 

argument  of  perigee  at  epoch 

Y 

(GAMMA) 

rate  parameter  of  nominal  anomalistic 
mean  motion  n 

T1 

(ETA) 

rate  parameter  of  the  major  axis  of  the 
nominal  ellipse 

T 

(TAU) 

rate  parameter  of  the  line  of  nodes  of 
the  nominal  plane 

no 

(CN) 

mean  motion  at  epoch 

n 

(CNBAR) 

nominal  anomalistic  mean  motion 

P 

(CP) 

semi-parameter  of  the  conic  at  epoch. 

d.  INPUTS 


ENTRY  I:  IENTRY  =  1  (initial  entry)  requires  t_q,  x_ q,  p,  1T>  IARG 


where 


IARG  =  1  time  is  independent  variable 

=  2  true  anomaly  is  independent  variable 
=  3  eccentric  anomaly  is  independent  variable 
ENTRY  2:  IENTRY  =  2  (normal  entry)  requires  and  A 
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jf.  •(; 

it  ; 

If  i 


il 


ir 
'  0 

**  \ 
}  ■ 


OfO  . 

(> 


<3 

!■' 


e. 


f. 


(100) 


OUTPUTS 

ENTRY  1:  IENTRY  =  1  gives  no  output  parameters 

ENTRY  2:  IENTRY  =  2  gives  nominal  trajectory  point  at  time  t  (r,  r,  r) 
COMPUTATIONAL  PROCEDURE  AND  EQUATIONS 
Test  IENTRY  =1 


Yes,  go  to  (100) 
No,  go  to  (200) 


Set  2ir 

=  2.  X  i r 

(SETUP  ENTRY) 

GMU 

=  H- 

Compute  initial  orbital  parameter  set 

ROSQ 

_  2 
-  r 

o 

=  r  •  r 
—  o  — o 

RO 

=  r 

o 

-  01/2 

RORDO 

=  r  r 
o  o 

~  x  •  r 
—  o  —  o 

VOSQMU 

-  *  2, 
-V1 

“  (i<,  '  41711 

RECIA 

=  <2/ro) 

•  2, 

~  ro  ^ 

CA 

=  a 

=  1 /RECIA 

CETERM 

=  C 

e 

=  1  -  (rQ/a) 

SETERM 

=  S 

e 

=  * 

0  0 

CE 

=  e 

=  (C  2  +  s  V^2 

e  3  1 

UOVEC 

=  U 
—  o 

-  r  /  r 
—  0  0 

WOVEC 

=  (£0xi0)/|r0: 

VOVEC 

=  V 
—  o 

=  V/  X  U 
—  o  —o 
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SINCI 

=  sin  i 

-  [u  2  +  v  2] 1/2 

oz  oz 

C0SCI 

=  cos  i 

=  f(u  +V  )2+(U  -V  )Z)lte 

ox  oy'  '  oy  or  1 

Cl 

=  i 

=  tan-*  (sin  i/cos  i) 

SINUO 

=  sin  u 

o 

=  U  /sin  i 
oz 

C0SUO 

=  cos  u 

o 

=  V  / sin  i 
oz 

cuo 

=  uo 

-  tan  1  (sin  uQ/cos  uQ) 

SINLO 

=  sin  i 

o 

=  (U  +  V  )/(l  +  cos  i) 

'  ox  oy 

C0SLO 

=  cos  i 

o 

=  (Uoy  -  Vox)/(1  +  COS  i} 

CLO 

=  l 

=  tan-  *  (sin  l  / cos  l  ) 

o 

o  o 

C0 

=  n 

=  i  -  u 
o  o 

smc0 

-  sin  J2 

C0SC0 

=  cos  & 

CVTERM 

=  Cv 

=  (p/rQ)  -  1 

SVTERM 

=  SV 

=  (p/p)1^2 

SINYO 

=  sin  v 

0 

=  Sy/e 

C0SVO 

=  cos  V 

o 

=  Cy/e 

YOARG 

=  V 

o 

=  tan-*  (sin  v  /cos  v  ) 
o  o 

ARGPER 

=  to 

II 

p 

o 

1 

o< 

Compute  initial  rate  parameter  set 

GAMMA  =  Y  =  -3/2  (a  /a J2  (a/r)3  (1-3  sin3  i  sin^  u  ) 

6  6  0  O  O  y 

ETA  =  t]  =  3/4  .T2  (ae/aQ)2  (1  -  e2}  “2  (4  -  5  sin2  iQ) 

TAU  —  t  ~  -3/2  J2  (ae/ao)2  (1  -  e2)  2  cos  iQ 
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i  0  0 


[Cl 

0  cos  i  sin  * 

_0  -sin  i  cos  i 

[a31 

ID]  [C]  [Aj]  [C“]  [D“] 

[A4] 

[d]  [cj  Ia2J  fc'J  [D“] 

[a51 

IAj]  IA3) 

fA6l 

lA2]  IA4] 

[A7] 

[a2]  Ia31 

[AgJ 

[Aj]  IA4] 

Exit 

(200)  Test  IARG  =  1,  2,  or  3 

If  equal  to  1,  go  to  (300) 

If  equal  to  2,  go  to  (400) 

If  equal  to  3,  go  to  (500) 

(300)  Time  is  the  independent  variable  (A  =  elapsed  time  from  epoch). 

Compute  mean  anomaly  difference  M  -  Mq 

DM  =  nA 
Call  KEQS 

Enter  with  DM,  C  ,  S  ,  e , 

6  6  K 

Exit  with  DE  =  E  -  Eq,  eccentric  anomaly  difference 
Compute  true  anomaly  difference  $  =  v  -  vq 


(301)  SIND.E  =  sin  (E  -  Eq) 

COSDE  =  cos  (E  -  Eq) 

R  =  r  =  a  1 1  -  C  cos  (E  -  E  )  +  S  sin  (E  -  E  )] 

L  e  o  e  o  J 

C0SPHE  =  cos  $  =  1  -  a2  (1  -  e2)  (1  -  cos  (E  -  Eq))  /  rrQ 

SINPHI  =  sin  $  =  [a2  (1  -  e2)1^2  /  rrQ] 

l(M  -  Mq)  -  (E  -  Eo)  +  sin  \E  -  Eq)] 

PHI  =  $  =  tan  *  (sin$  /cos  $) 

The  monotonicity  adjustment  for  $  »  f  -  f  is 

If  $  >0,  then  $  =  $  +  AE  -  MOD  (AE,  2 it),  or 
K  |$|  -  |ae|  >  tr»  then  $  =  $  -  2  tt 
If  $  <  0  (AM  <  0),  then  $  =  $  +  AE  -  MOD  (AE,  -2ir),  or 
If  j$j  -  |ae|  > ttj  then$  =  $+2ir 

Go  to  (600),  ,, 

(400)  True  anomaly  is  the  independent  variable  (A=  4>  ~  v  -  vQ). 

Compute  eccentric  difference  E  -  Eq 


R  =  r 

=  p  [  1  +  C  cos  $  -  S  sin  $]" 1 

r  *•  \r  V 

SINDE  =  sin  (E  -  Eq) 

1/2 

=  r(ap)  '  sin$  -  r/p  (1  -  cos$)  Sg 

C0SDE  =  cos  (E  -  E  ) 

0 

-  1  -  rrQ  (ap)-1  (1  -  cos  4>) 

DE  =  (E  -  Eq) 

=  tan^^sin  (E  -  Eo)/cos  (E  -  Eq)) 

Compute  me?ji  anomaly  difference  M  -  M 
*  ■  0 

DM  =  (M  -  M  ) 

o 

-  (E  -  E  )  +  2  S  sin2  (E  -  E  /2) 
o  e  o 

. 

-  C  cos  (E  -  E  ) 
e  .  o' 
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Compute  corresponding  time  parameter  (t  -  tQ) 

DT  =  (t  •  t  )  =  (M  -  M  )  /  n 
o  o 

Go  to  (600). 

(500)  Eccentric  anomaly  is  the  independent  variable  (A=  E  -  Eq). 

Compute  mean  anomaly  difference  M  -  M0 

DM  =  (M  -  Mq)  =  (E  -  Eq)  +  2  Sg  sin2  (E  -  Eq/2)  -  Cfi  cos  (E  -  Eq) 
Compute  corresponding  time  parameter  (t  -  t^) 

DT  =  (t  -  tQ)  (M  -  Mq)  /  n 
Go  to  (301). 

(600)  Compute  reference  trajectory  point  <on  the  conic  at  time  t 

(r>  £>  £).  Evaluate  the  F  and  G  function  coefficients 

F  =  Fj  r  ccs  $  +  F.,  r  sin  $ 

G  =  Gj  r  sin$ 


Evaluate  the  scalar  coefficients  a.  (j=l,  ....  6 )  of  the  rotation 
matrix  [A] 

a^  =  sin  (t$) 

a2  =  -2  [sin  (t$/2)]2 

a^  =  sin  (n$) 

a4  =  -2  [sin  (*l$/2)]2 

> 

a5  "  aia3 
a6  =  a2  a4 
a7  a2  a3 
a8  =  al  a4 


(B-8) 
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Evaluate  the  rotation  matrix  [a] 


[A] 


=  £, 
k  J 


(«>  L-A..3 


Evaluate  the  [F]  and  [G]  matrices 

[F]  =  f[i  +  a] 

[Gj  =  G  [i  +  A  3 

(620)  Evaluate  the  reference  position  vector  £ 

c  r  =  [Fj  r  +  [G]t 
—  — o  — o 

Evaluate  the  F  and  G  function  coefficients 
AE  =  na/r 

r  =  a  AE  [C  sin  (E  -  E  )  +  S  cos  (E  -  E  )] 


•sinti  =  -  a2  (1  -  e2)jain  (E  -  EQ)  AErrQ  -  y  (1  -  cos  (E  -  Eq»] 

<rr/ 


COS  $$  = 


a2  (1  ^  e2)l/2" 
— ; — 72 - 

(n  -  AE  (l  -  cos  (E  -  E  ))] 

-  <«o> 

1  \  o  •  i 

>1  vl 


-  [(M  -  M  >  -  (E  -  Eo)  +  sin  (E  -  E( 

$  =  [(-  sin$6)2  +  (cos 

F  =  Fj  r  cos  $  +  Fj  r  (-  sin  $$)  +  F^  r  sin  $  +  r  (cos  $6) 

G  =  Gj  r  sin®  +  Gj  r  (cos$$) 

Evaluate  the  scalar  coefficients  a.  (j  =  l»  ...»  8)  of  the  rotation  matrix 

[a] 

a.,  =  T®  (a_  +,  11 

i  *  c 

^2  =  al^ 


-86- 


a3  =  -n  6  (a4  +  1) 

a4  =  ii#  {-  a3) 

a5  =  aia3  +  a2a3 

a6  =  a2  a4  +  a2  a4 

a7  =  a2  a3  *  a2  a3 

a8  =  al  a4  +  al  a4 

Evaluate  the  rotation  matrix  Ia] 

8 

1X1  "X  (9,  4)  lAj] 

i=i 

Evaluate  the  [fJ  and  [Gj  matrices 

If]  =  [f[i  +  A]  +f[a] 

16 J  =  [g[i  +  A]  +  G  [A] 

(660)  Evaluate  the  reference  position  vector  £ 

r  =  IF]  r  +  16]  Jr 

—  — o  — o 

Evaluate  the  different5, al  operator  (  ;  ) 

(r  cos#  ;  a.)  (j  =  1,  . . . ,  8) 

J 

(r  sin#;  a.)  (j  s  1,  . .  8), 

J 

where  (p  ;  q)  =  2  (dp/dt)  (dq/dt)  +  p  d2q/dt2 
Evaluate  the  [F]  and  16  ]  matrices 


F 


(r  cos  $;a.)  +  F2  (r  sin  $;a.)][A.}=  {F;[A]) 
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that  is, 


AE  =  -  f(E  -  Eq)^/  f'(E  -  Eq)^  (Newton-Raphson  method) 


Use  the  first  approximation  of  (E  -  Eq)  given  above  as  an  initial  guess, 

where 

f  (E  -  Eq).  =  (M  -  Mq)  -  (E  -  Eq).  -  2  Se  sin2  [  (E  -  Eq)/2] 

+  C  sin  (E  -  E  ). 
e  o'i 


and 


give 


f  (E  -  E  )  -  -1  -  2  S  sin  [  (E  -  E  )  /2] .  cos  [(E  -  E  )/?,]. 

o  e  g  i  o  i 

+  C  cos  (E  -  E  ). 
e  oi 

=  -1  -  S  sin  {E  -  E  ).  +  C  cos  (E  -  E  ). 

e  o  i  6  o  i 

AE  _  (M  •  M  )  -  (E  -  E  ).  -  2  S  sin2  [(E  -  E  J/2J.  +  C  sin  (E  -  E J 

1  +  S  sin  (E  -  E  ).  -  C  cos  (E  -  E  ). 
e  o'i  e  o  i 


Until  AE  <  € 


(10”8), 


this  procedure  yields  a  solution  (E  -  E  ). 
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